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Abstract. We present in this paper a simplification of the dune model proposed by Sauermann et al. which 
keeps the basic mechanisms but allows analytical and parametric studies. Two kinds of purely propagative 
two dimensional solutions are exhibited: dunes and domes, which, by contrast to the former, do not show 
avalanche slip face. Their shape and velocity can be predicted as a function of their size. We recover in 
particular that dune profiles are not scale invariant (small dunes are flatter than the large ones) , and that 
the inverse of the velocity grows almost linearly with the dune size. We furthermore get the existence of a 
critical mass below which no stable dune exists. However, the linear stability analysis of a flat sand sheet 
shows that it is unstable at large wavelengths and suggests a mechanism of dune initiation. 

PACS. 45.70.-n Granular systems - 47.54.+r Pattern selection; pattern formation 



1 Introduction 

' The beauty of the crescentic barchan dunes have recently 
[ attracted the interest of physicists for a better understand- 
• ing and modelling of sand transport, as well as ripples 
' and dunes formation and propagation. E. Guyon had this 
, witty remark about them: 'barchans are our drosophila', 
■ which means that beyond the scientific and fundamen- 
tal works on these dunes, we all keep in mind that such 
1 studies may lead to potential applications in the fight of 
saharan countries against sand invasion. One of the first 
reference work in the field is certainly the famous book 
of Bagnold which dates back from 1941 [|lj. Since then, a 
great effort of measurement and modelling has been done 
which we have reviewed in details in the first part of these 
twin papers. 

Our aim here is to discuss and model the selection of 
two-dimensional dune shape and velocity. For that pur- 
pose, we will simplify the model proposed by Sauermann 
et al. |§,||,D . We will show that although rather severe ap- 
proximations, we are able to recover their main results, 
in particular that dune profiles are not scale invariant, 
and that the inverse of the velocity grows almost linearly 
with the dune size. Besides, analytical expressions of dome 
and dune propagative profiles can be obtained, but whose 
coefficients have to be numerically computed. We further- 
more get the existence of a critical mass below which no 
stable dune exists. The apparition of dunes can however 
be understood with the linear stability analysis of a flat 
sand sheet which is unconditionally unstable towards large 
wavelengths perturbations. 



The paper is organized as follows. Section g is devoted 
to the equations of the model. The linear stability of a uni- 
form sand sheet is treated in section]^. In section ^ we sim- 
plify further the equations and show what is the general 
shape of the purely propagative solutions of the model. 
The specific case of domes and 'actual' dunes are discussed 
in sections ^ and |^ respectively. At last, we conclude with 
a discussion of the relevance of these results, and the possi- 
ble extension of the model to three-dimensional situations 
and dynamical studies. 

2 Basic equations 

We wish to give a description of the shape and evolution 
of two-dimensional dunes in terms of two fields: the profile 
h{x,t) and the volumic sand flux q{x,t) which is the vol- 
ume of sand transported through an inflnite vertical line 
per unit time, x denotes the horizontal coordinate and t 
is the time. We are going to write a set of three equations 
for these quantities in order to include into the model 
(i) the mass conservation, (ii) the progressive saturation 
of sand transport and (iii) the feedback of the topogra- 
phy on the sand erosion/deposition processes. Although 
we shall keep in this paper to two-dimensional situations 
which correspond to transverse dunes (invariant in the di- 
rection perpendicular to the wind), our ultimate goal is of 
course to be able to describe three dimensional dunes and 
barchans in particular. 

Because barchans in dune flelds - so-called 'ergs' - 
organize themselves like gooses or ducks during their mi- 
gration flights, we named this class of models for purely 
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Fig. 1. When the wind is blowing on a patch of sand, the flux 
of transported sand get saturated after a typical length Isat 
which is almost independent of the wind shear velocity u* - 
see part 1. 



graphical reasons the modellings. This denomination 
includes the approach of Sauermann et al. as well as the 
different variations and simplifications we derived in this 
paper from their work. At this stage of the modelling, we 
are however far from being able to take into account such 
'interactions' between dunes which are necessary to ex- 
plain the Cq spatial organisation, and we shall focus on 
isolated objects only. 

A simple balance calculation shows that the erosion 
rate —dth is directly related to the divergence of the flux 
q. This gives the common continuity equation: 



dth + dxq = 0. 



The saturation effect of the sand transport has been 
already evoked in the first part of these twin papers: con- 
sider a patch of sand on which the wind is blowing - see 
figure |l|. The flux of transported sand q first increases and, 
because of the feedback of the grains on the wind velocity 
profile, get saturated after a typical length Isat ■ In the first 
part of the paper, we showed that I sat = £,d Psand/ Pair {d 
is the grain diameter, p the densities and ^ a non dimen- 
sional prefactor), i.e. is almost independent of the wind 
shear velocity - slow logarithmic dependencies hidden 
in ^ only. This phenomenon has been reported and stud- 
ied by several authors, e.g. [0J|]. The real shape of q{x) is 
certainly more complicated than the one drawn on figure 
^. In particular, oscillating or overshooting features were 
reported in |^ when q reaches its asymptote. However, 
what is important for our purpose is only that a satu- 
rated value qsat is reached after a length Isat- This space 
lag is satisfactorily described by the following equation: 



Qsat 



^ sat 



(2) 



This charge equation can be seen as a simplification 
of that proposed by Sauermann et al. in their continuum 
saltation model ||^ . An important remark is that this equa- 
tion is valid only if some grains are available on the sand 
bed. On a firm soil indeed, the fiux cannot increase to be- 
come saturated. As suggested by Peer and Hakim the 
right hand side of the relation (0) must be therefore mul- 
tiplied by some matching function which quickly tends to 
zero when h is decreased below, say, d the grain diame- 
ter, and which is equal to unity above this value - the 



altitude of the firm soil is /i = 0. Then, the equation (g) 
becomes non linear but no boundary conditions have to 
be specified at the edges of the sand covered region. For 
example, if the matching function tends to zero like h, the 
dune will always keep a thin sand sheet at its back. But if 
it varies as ^/h, the dune will have a finite extension and 
will join the firm soil with an horizontal tangent. We shall 
ignore at present these subtleties but keep them in mind 
to invoke them later when necessary. 

Another important remark is that the time scale on 
which the dune profile h evolves is incomparably larger 
than that of the sand flux q. We then assume that q adapts 
its profile instantaneously according to equation (H) and 
makes h change slowly through equation (|l|). Therefore 
any term dtq is irrelevant in this modelling. 

The saturated fiux qsat is uniform for a flat sand bed 
only. To the first order, the saturated flux qsat is a function 
of the local shear stress t = PairU^ which itself depends 
- non locally - on the topography: basically, bumps and 
upwind slopes get more eroded than dips and downwind 
faces. A classical relationship between the saturated flux 
and the shear velocity that can be recovered with the scal- 
ing arguments of the part 1 of the paper is: 



Qsat oc 



id 9 



(3) 



In principle, such a relationship is valid far from the veloc- 
ity threshold Uthr under which no sand can be eroded by 
the wind, i.e. qsat — for < Uthr- Refined formulas can 
be obtained which essentially smooth the step from to 
the previous asymptotic expression such as that obtained 
in part 1: 

Pair "U* / 2 2 \ f a\ 

qsat OC [U^ - Uthr) ■ 



Psand Q 

In the whole range —uthr < m* < Uthr, Qsat is null so that 
in practice, qsat cannot become negative on a dune. This 
condition will be used in section |[ 

To close the equations, we have to explicit the spa- 
tial variations of the turbulent wind velocity due to the 
dune profile. The simplest model which verifies the ba- 
sic requirements (see part 1) is certainly the perturbative 
calculation by Jackson and Hunt Neglecting loga- 

rithmic scale dependencies, Kroy et al. ^] have extracted 
the main features of their work by expressing the shear 
velocity as: 



= 1 



A [ — dxh{x 

J TTS 



s) + Bdxh{x), (5) 



where J7* is the shear velocity exerted on a flat bed. More 
precisely, as discussed in details in the part 1, A and B in 
principle depend on the size of the dune D with IhD/zq 
factors, where zq is the roughness of the sand surface. For 
D varying between 20 and 200 m and a roughness of order 
of the grain size, such a logaritmic factor does not change 
by more than 20% over the whole range. This justifles the 
fact that it is reasonable to take constant effective values 
for the coefficients A and B. Several further important re- 
marks must be made on this expression. First, it must be 
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noted that the convolution integral acts on dxh roughly 
like a derivative, leading to a term which encodes curva- 
ture effects. But this curvature is dimensionless and thus 
does not depend on the dune size - in other words, this 
term can be seen as a curvature rescaled by the dune size. 
It reflects the observation that the wind velocity increases 
on bumps (negative curvature) and decreases on hollows 
(positive curvature). Second, it is a non local term, mean- 
ing that the shear velocity depends on the whole shape 
of the dune. Of course sharp variations of the dune profile 
will also have a strong local effect. At last, the second term 
simply takes into account slope effects: positive slopes are 
more eroded than negative ones. Again, this term does not 
introduce any new lengthscale. 

Expression (^) can be used to close up the set of equa- 
tions as was done by Kroy et al. in Q] who used besides 
a more sophisticated - and non-linear - charge equation 
than (^) . It is useful to simplify further the equation link- 
ing Qsat to h without loosing too much physics, in order to 
let more analytical developments. The linear expansion of 
the expression (^) rewritten in terms of qsat oc gives: 

^f^^l + ^A f-d,h{x-s) + lBdxh{x), (6) 

Qsat 2 J TTS 2 

where Qsat = qsatiU^:) is the saturated flux on a flat sand 
bed submitted to a shear velocity C/,. 

Furthermore, we shall use the saturation length l^at 
and flux Qsat to make our variables dimensionless. Thus, 
for a given wind shear velocity, all relevant scales of the 
problem are fixed. For instance, Qsat /I sat is the veloc- 
ity scale, llg^^/Qsat the time scale, Qsat/l'^at ^he frequency 
scale, etc. Note that the strength of the wind is completely 
encoded in these rescalings. 

3 Stability of a flat sand bed 

Before going further in the modelling of dune shape and 
propagation, we wish to investigate the problem of dune 
initiation. Indeed, there are two striking field observations. 
First, no persistent barchan dunes exist smaller than say, 
1 m high, 20 m large and 20 m long. Second, any small 
conical sandpile blown by the wind disappears even when 
a sand supply is provided. Then, how can barchan dunes 
appear? 

To investigate this problem, we have integrated numer- 
ically equations (|l]j2[||) with A = 6 and i? = 4. Two initial 
conditions were tested. First, a small triangular sandpile 
at the repose angle is prepared on the firm soil (figure ^ . 
It can be seen that it is rapidly eroded and disappears, 
as observed in the field and in wind tunnel experiments 
- see part 1. Second, we look at the evolution of a thin 
sand sheet (figure |3|) disturbed by a flat bump. This initial 
conditions mimics a sand beach on which sand is deposed 
by water. It can be seen that the bump propagates down- 
wind and induces a strong erosion of the sand bed in front 
of it. A second bump nucleates from the initial pertur- 
bation which itself induces a strong erosion in front of it, 
and so on. After some time, a series of growing oscillations 



li(x,t) 
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Fig. 2. Numerical integration of equations (|ljy|6|) for the evo- 
lution of a small triangular sandpile on the firm soil, initially 
1.3 Isat high and at the repose angle. h{x, t) and x are in units of 
laat. The time between two profiles is 0.1 in units of Isat /Qsat- 
For legibility, the profiles are translated vertically from time 
to time. The grey filled region shows the available amount of 
sand at t = 0. The field observation that such a small sandpile 
disappears when blown by the wind is then recovered in the 
model. 



is generated. The amplification of this phenomenon stops 
when the oscillations eventually reach the firm soil from 
which no sand can be eroded, and/or when recirculation 
bubbles - see section ^ - appear and make these bumps 
interact. As a conclusion, depending on their spatial ex- 
tension, small sand bumps either disappear or grow and 
initiate dunes. 

It is then instructive to make the linear stability anal- 
ysis of a flat sand bed. Let us consider an infinite uniform 
sand bed blown by a uniform wind. The sand flux is every- 
where saturated: g = 1 (in units of Qsat)- To investigate 
its stability, we can consider, without loss of generality a 
small perturbation of the form: 

/i(a;,i) = i?e'"*-'"*+**-'^, (7) 
g(a;,i) = l-^Qe'"*-"^*+^'=^ (8) 

where Q and H are related one to the other by the con- 
servation of matter, 

{iu-a)H = ikQ. (9) 

^From the relation (||), we get the following expression 
for the saturated flux: qsat = 1 + 3/2(|fc|A -|- ikB)h. Once 
replaced in the saturation equation, it gives: 

(1 + ikhat)Q = \{A\k%at + Blk)H. (10) 
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Fig. 3. Numerical integration of equations (|1[|2|J6|) for the evo- 
lution of a flat sand sheet of width Isat over the firm soil dis- 
turbed by a small hump. h{x, t) and x are in units of hat- The 
time between two profiles is 0.5 in units of Isat /Q sat - For legibil- 
ity, the profiles are translated vertically from time to time. The 
grey filled region shows the available amount of sand at t = 
and at the final time. The perturbation propagates downwind 
and growing oscillations nucleate from it. The hollow created 
in front of the initial bump soon reaches the firm soil. 



Combining equations (|9|) and ( |10D we finally obtain 

_ 3fc2(B- A|fc|) 
~ 2(l + /c2) ' 

_ 3k\k\{A + B\k\) 
" 2(l + /c2) ■ 



(11) 
(12) 



Note that a more complicated relation between qsat and 
would only affect the prefactor, i.e. the time scale, in 
this calculation. 

The growth rate a, shown for A — 6 and _B — 4 on 
figure 1^ (top), is positive for small wavenumbers (fc < 
B/A) and negative for large wavenumbers (fc > B/A). A 
flat sand bed thus exhibits a large wavelength instability 
which can explain the initiation of dunes. But it is stable 
towards small wavelength disturbances, so that a small 
sandpile on a firm soil is quickly eroded. 

The fastest growing mode - that which maximises a 
- is obtained for 3k + = 2B/A. Neglecting the k^ 
term, we obtain a still good approximation of the most 
unstable wavelength: A = 27T/k ~ inlsatA/ B. For A = Q 
and i? = 4 we get A ~ lAlsat- For I sat — 9 m we get a 
wavelength which is reasonable compared to what is ob- 
served on transverse dune fields in deserts. Under water 
Isat rather scales on the grain size and we have Isat — 1 cto 
which gives again a good estimation of what has been 
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Fig. 4. Results of the linear stability of a uniform sand bed 
blown by the wind. Top: growth rate a, rescaled by Qsat/lsat, 
as a function of the disturbance wave number k, rescaled by 
i/hat- The sand bed exhibits a large wavelength instability. 
Bottom: group velocity c of disturbances, rescaled by Qsat /hat, 
as a function of the disturbance wave number k, rescaled by 
i/lsat- For small wave numbers, the velocity c increases linearly 
with k and thus decreases as the inverse of the wavelength. 



measured by Betat et al. |^ or Andersen et al. jl^ for 
example. 

Finally, it can be seen from equation (^2|) that distur- 
bances propagate downwind. We can compute the group 
velocity of these surface waves: 



duj _ 2,[A\k\+ Bk^ii + k"^)] 
'dk ~ 2(l + fc2)2 ■ 



(13) 



c is plotted on figure ^ (bottom) for A = 6 and B — A. 
For small wavenumbers fc, the group velocity c increases 
linearly with k. This means that for asymptotically large 
wavelengths A, the propagation speed scales as Qsat/^- We 
thus recover in the limit of large bumps, the scaling pro- 
posed by Bagnold. Wc also see that the velocity deviates 
significantly from this law, when the wavelength becomes 
comparable to the saturation length Isat- This is also the 
case for actual dunes - see part 1. 



4 A simple 2d modelling 
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4.1 Simplified equations 

The expression (|^) can be kept as it is, mixed with the 
mass conservation and the charge (||) equations, and 
treated in Fourier space, but the idea is rather to replace 
the convolution expression by a simpler anzats. Since it 
is a 'scaleless and non-local curvature', we can replace it 
by —Ddxxh, where D is the dune length. D will be deter- 
mined by the boundary conditions. We then end up with 
the following set of - linear - equations of our model: 

dth + d^q = 0, (14) 
dxq = qsat - q, (15) 
qsat = 1 - aDdxxh -I- I3dxh. (16) 

a and (3 are two of the parameters of the model - a third 
one, /if,, will be introduced later on. They could be in 
principle measured independently. The a term has a sta- 
bilizing role on small perturbations, while the (3 one make 
them growing. The rescalings h <— ah and t <~ at (with 
q and x unchanged) shows that only the ratio /3/a is ac- 
tually important. As stated in the figure captions, most 
of the curves of the paper have been plotted with a — 1. 
and /3 = 4. At last, as we already mentioned, the charge 
equation (|l^) is valid only if there is some sand to be 
eroded h > 0. Note besides that, to be fully consistent, 
qsat should remain positive everywhere. 

4.2 Boundary conditions 

Before getting into the process of solving this system of 
equations, we need to specify what the boundary condi- 
tions at the edges of the dune are. This rises several im- 
portant questions. The first thing to notice is that, when 
we transformed the relation (|^) which links q^at to the 
profile h with a convolution term into the simpler expres- 
sion ([l6|), we have incremented the order of the differen- 
tial equation. Therefore, such a simplification requires an 
additional boundary condition than what is needed to in- 
tegrate equations (^,|l^,^5|) for example. First order equa- 
tions like the later need to know the profile h at some posi- 
tion. However, when using the second order equations (|lj- 
[l6| ), the slope h' should be also specified. As we already 
mentioned, one way to solve this problem is to regularise 
equation ( p^ by a non linear factor when h tends to zero 
which gives conditions on h(x) and h'{x) as x — s- —oo, see 
discussion after equation (H) in section ^. 

An alternative is to consider the surface profile - the 
firm soil as well as the dune - as a whole. This is pre- 
cisely what is done when the convolution (|^) is used. If 
the slope h' was not continuous at the dune edge, in some 
sense, h" would be infinite. Then, from equation (p^), the 
saturated flux would formally tend to — oo over a region 
around this edge, i.e. would be forced to be null because it 
cannot be negative. Thus in this region, only sand deposi- 
tion is permitted. This means that a discontinuity in the 
slope at the dune edge, immediately reacts to prevent the 
motion of this boundary point, which starts moving again 
only when the profile becomes flat, as shown on figure ^. 
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Fig. 5. The numerical integration of equations for an 

initial sand pile making a finite angle with the firm soil shows 
that the profile of its upwind left edge quickly evolves to get a 
horizontal slope which smoothly continue the soil. h{x, t) and 
X are in units of l^at- The displayed profiles corresponds from 
left to right to the times t = 0, 0.25, 0.5, 0.75 and 1.5 in units 

of hat/Qsat- 



Therefore, considering that the wind and thus the satu- 
rated flux are not sensitive to the transition from the firm 
soil to the dune back, gives a 'natural' upwind boundary 
condition for propagative solutions: the profile h as well as 
its slope h' should vanish at the upwind boundary of the 
dune. This argument actually also applies at the down- 
wind edge, and we shall use it in the dome and dune next 
sections. 

In the sequel, we are going to look for propagative so- 
lutions, i.e. functions of the type h{x — ct) and q{x — ct). 
We shall see that the explicit form of these functions can 
be found for given upwind conditions. However, such solu- 
tions are parametrized by coefficients such as this length 
D or the propagative velocity c which must be fixed by 
the right - downwind - conditions. Two kinds of right 
boundary conditions will be considered, leading to so- 
called 'dome dunes', i.e. without avalanche shp face, or 
'actual' ones for which a 'recirculation bubble' will be in- 
troduced. 



4.3 General form of the propagative solutions 

For functions of a:; — ci, the continuity equation can 
be easily integrated and gives 

g = go + ch. (17) 

We now describe everything in the propagating frame ref- 
erential, and rename x as the new space coordinate. We 
look for isolated propagative objects. The point a; = is 
the beginning of the dune where h = Q. go is thus the in- 
coming sand flux - the sand supply. In the region x < 0, 
no grains are available on the ground: h = Q and q — qo 
everywhere. In the region a; > 0, using the relation be- 
tween g and h ( pT[ ) in the charge equation ( p^ with q^at 
given by its expression ( |T6| ) , we get an ordinary differential 
equation for the dune profile h: 

l-aDh" + {(3-c)h' ~ch-qo = 0, (18) 
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Fig. 6. Longitudinal profile of a barchan dune at Negrita 
beach, southern Morocco. Black and white dots correspond 
to independent measurements of the same dune. The solid line 
correspond to the best fit by a function of the form ([f9|). 



where h' and h" denote the first and second derivatives of 
the dune profile. 

Under the two conditions /i(0) = and h'ifi) — 0, the 
solution of equation (nsh is 



h{x) 



1 - go 



+ {^^ sin{kx) - cos{kx)^ e'"' , (19) 



where the coefficients s and k are given by 



k = 



2aD ' 
1 



2aD 



^/AcaD -{[3- cf 



(20) 
(21) 



Figure ^ shows the comparison between the central 
longitudinal profile of an actual barchan and the theoret- 
ical form (n9|) which is parametrised by fc, s and c. It can 
be seen that the shape of the dune as well as the bound- 
ary conditions are well captured by the model. In fact, an 
oscillatory function as (|l9|) does not look like an isolated 
dune yet. In particular, this h{x) takes negative values. 
Besides, as already mentioned, the two coefficients D and 
c are not determined yet. In the next sections we show 
how it is possible to cut this solution at some point with 
an adequate right boundary condition in order to get a 
genuine propagative profile. 



5 Domes 

In this section, we look for so-called 'dome' solutions, i.e. 
dunes which do not show any avalanche plane. This is 
possible if the local slope of the dome is everywhere not 
steeper than some threshold /Ltf,. In the next section, it 
will become clearer that this threshold corresponds to the 
slope at which the wind stream lines detach the dune pro- 
file, creating a backward wind flow - or a 'recirculation 
bubble' - which leads to a slip face. As will be explained 
in the conclusion, these dome solutions may play an im- 
portant role when extending the present model to three 
dimensional situations. One of the major results of this 
section is that the shape of the dome, when it exists, is 
selected by the value of the incident flux q^. Because we 
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Fig. 7. Shape of the dome solution for different values of the 
sand flux go- The data have been computed with a = 1. and 
13 = A. Note that these domes are actually very flat: for clarity 
the scale of the vertical axis is much larger than that of the 
horizontal one. 



look for purely propagative solutions, go of course is also 
the flux of the sand which leaves the dune. 

Let us first take a very simple example. For c — [3 the 
coefficient s defined in equation (|2^) vanishes, and the 
solution (|l9| ) reduces to h{x) = [1 — cos(fca;)] (1 — qo)/f3, 
where k = [/3/(aZ))]-^/^. As the general solution ( p^ ) was 
constructed for, h and h' vanish at a; = 0, but also at a; = 
211 /k. This value sets the length of the dune: D = 47r^a//3. 
Interestingly, this is precisely the cut-off Ac of the linear 
stability analysis. As discussed in the previous section, this 
simple solution satisfies the 'natural' downwind boundary 
conditions h{D) — and h'{D) — 0. Another important 
restriction is that on the saturated flux qsat which must 
remain positive everywhere. As a matter of fact, there is 
a minimum value of 50 below which qsat takes negative 
values and makes this solution inconsistent. It is easy to 
show that this lower bound for go reads 



1 



(22) 



For the values a = 1. and /3 — 4., we get q^ ~ 0.16. Two 
such dome proflles for go = 0.5 and go = 9o displayed 
on figure 0. At last, note that the steepest slope of this 
solution is h'^^^^ = For go = ^"^^ ^^'^ same 

numerical values as above, it gives |^mj„| — 0.13. For any 
slope threshold fih larger than |/imj„|, the whole set of 
these solutions such that go < go < 1 is acceptable. For /if, 
below this (actually rather small) value however, go must 
be larger than gf; = 1. — 27ra/Zb, and no fully consistent 
dome solution can be constructed for a smaller incident 
sand flux. In the sequel, the typical value of /i^ that we 
shall use will be fii, = 0.25. 

How can we construct a dome solution for go < go? 
Suppose the values of the velocity c and the length D 
are given, from equation ( |l9|) we can compute h and its 
derivatives as well as qsat- Because go is smaller than gp, 
qsat will reach zero at some position x = L < D. Negative 
values of qsat £^re not permitted, and we therefore set it 
to zero for x > L. Then, equation (^|) is very easy to 
integrate and gives and exponential branch for the flux 
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Fig. 8. Definition of the parameters of the dome, a; = is the 
reference point of the dune. H is its height. For qo < Qq, its 
length D must be cut into three regions of size L, R and Y. 
As explained in the text, in the middle region (dahed line) the 
saturated flux vanishes and the profile is a branch of exponen- 
tial which matches with the two parts of the full solution. Hr 
and Hy are the heights at the two sides of this central region. 



q{x), which, using the linear relation (|T^ between q and 
h leads to: 



hr{x) 



(23) 



The subscript r is used to avoid any confusion with h{x) 
given by equation (p^, but both are part of the same 
dome solution. In particular, the dome must be smooth 
and in this expression Hr — h{L) - we shall come back in 
the next paragraph on the continuity conditions at x = L. 
We cannot end the dome solution with hr for two reasons. 
First, we argued in the previous section that the dome 
profile must ends with a horizontal slope, which is not 
possible with an exponential function at finite distance 
D. Second, if one compute what would be the saturated 
flux gLt calculated from the profile hr with a relation 
like dlQ), one sees that it is negative as it should close to 
X = L(remember that qsat is set to zero) , but crosses zero 
at some other position x ^ L + R < D. From that point 
and for larger x we thus need to come back to the original 
profile h. Then, a natural way to end up the dome profile 
is to use 

hy{x) = h{D - x), (24) 

which is consistent if qsat{—Y) — 0, with L + R + Y — D. 
This choice ensures that both h{D) and h'{D) vanish as 
required. These three regions of the dome solution are 
illustrated on figure |^. 

Let us now be more explicit about the continuity con- 
ditions at the two matching points. At x = L, by con- 
struction of the relation (p^), the dome profile, the sand 
flux, and the saturated flux are continuous. Because equa- 
tion ([l5| ) holds everywhere, it implies that dxq-, and there- 
fore the slope h' are also continuous. By contrast, the 
curvatures h"{L) and hr{L) are different, and therefore 
llati^) 7^ qsat{L) = 0. The position x — L + R is defined 
by q^ati^ + R) = qsati—Y) = 0. At this point of course 
we do not want any step in the dome profile, such that 
Hy = hy{L + R) and hr{L + R) must be equal. Again, 
because of equation ( p^ ) and the fact that qsat has been 
built to be continuous at x = L + R, the continuity of the 
profile makes the slope continuous too. However this posi- 
tion is also the point where the pseudo saturated flux q^^^ 




6 9 
horizontal position x 

Fig. 9. Dome profile and its derivatives, as well as the satu- 
rated sand flux for go = 0.09. With this value of the incident 
sand flux, the dome steepest slope just reaches its minimum 
permitted value —fib at x = L. qsat is strictly zero between 
X — L and x = L -\- R and h is a branch of exponential. At 
X = L the curvature h" shows a discontinuity, but all quantities 
are continuous at x = L + R. 



crosses zero, such that the curvature of the dome profile 
is also continuous. All these continuity conditions can be 
shown on figure ^. 

To sum up, we have in practice four coefficients to 
determine: c, L, R and Y {D = L + R + Y) with the four 
non-linear following equations: 



qsat{L) = 0, 

qsat{-y) = 0, 
hr{L + R) = h{^Y), 
gL,(L + i?)=0, 



(25) 
(26) 
(27) 
(28) 



which can be done numerically. One can find such a solu- 
tion for any qa < q^. However, to be consistent with the 
next section, and as already mentioned in the c = (3 case, 
such a dome is acceptable only if its steepest slope is larger 
than —fib- This fixes a lower bound q^ under which there 
are no solution. With a = 1., /3 = 4. and fMb = 0.25, we 
find that <?§ ~ 0.09. The solution for this particular value 
of the incident sand flux is plotted on flgure |^, and the 
dome profile for go = 0.1 is compared, on figure |^, with 
two cases of the c — P solutions. 

When we have all coefficients of a particular solution, 
we can compute the mass of the corresponding dome: 



M = / dxh{x). 
Jo 



(29) 



For the c = (3 solutions, this mass can be simply expressed 
as M = 47r^a(l — qo)/ P'^. The interesting point is that the 
function A/(qo) decreases with its argument, see figure |o[ 
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Fig. 10. The main graph shows the mass of a dome selected by 
an incident flux qo. This function is a straight line for go > go 
which corresponds to the c = /3 solutions. The solutions for 
go < go are cut off at go — go where the slope at a; = L is equal 
to the threshold —ni, (big dot). The important point is that 
this curve M(go) is a decreasing function, which make these 
domes unstable to a change of go. In the inset, the inverse of 
the velocity 1/c is plotted against a typical size M^/^ of the 
dome. It is strictly constant for M < A/(go) and then almost 
straight up to the cut-off value Mc- 



* L R * 

Fig. 11. Definition of the dune and recirculation bubble pa- 
rameters, x = is the reference point of the dune. The brink 
- where the slip face and the recirculation bubble start - is 
located at x = L. The separaration streamline extends the 
dune profile, starting with the same slope p and curvature C . 
It smoothly reaches the ground a.t x = L -\- R. We note H the 
height of the highest point of the dune (crest). B is the dune 
height at the brink. H and B may coincide or not. We call ip 
the angle of the slip face. For the numerics, we use = 30°. 



Then, take a propagating dome under a flux go, and sup- 
pose that, for some reasons, the incident flux is suddenly 
increased by some small amount 5qo. The shape of the 
profile is such that the flux of the sand which leaves the 
dome is exactly go- Then, its mass gets larger. But a larger 
mass means that the flux of sand which can leave the dome 
is even smaller than before, which makes the mass of the 
dome increase again, and so on. In the same way, if the ini- 
tial flux 50 is suddenly decreased, the dome will loose more 



and more sand because less mass means a higher leav- 
ing sand flux. This negative feedback mechanism suggests 
that domes are unstable objects that may either quickly 
disappear or reach the point where their slope is steep 
enough to generate a recirculation bubble and create an 
avalanche slip face to become an actual dune. One could 
however observe them, for example under water with peri- 
odic boundary conditions which can stabilize such an 
instability. Then an interesting prediction on the velocity 
of these domes which could be compared to experiments 
is shown in the inset of figure ^ it is constant for small 
domes and 1/c grows almost linearly with their size M^/^ 
for larger ones. The dome length D behaves also very much 
like 1/c as a function of M^/^. 

6 Dunes 

It is possible to get another type of propagative profiles - 
transverse dunes - with a quite more sophisticated down- 
wind boundary condition, namely a recirculation bubble. 
It is a common field observation (see part 1) that the wind 
streamlines on a dune follow exactly the shape of its back 
profile but separate at the point where the avalanche slip 
face begins - the brink - and reattach further downwind, 
as shown on figure |ll|. This phenomenon creates an eddy 
in the 'shadow' of the dune, where the wind is much less 
strong than anywhere else. As a consequence, all the sand 
eroded on the back is deposited around the top of the slip 
face which avalanches when the slope becomes too steep. 
As explained in the first part of these twin papers, it is 
fortunate that an accurate description of these avalanches 
is not necessary due to the fact that they do not have any 
feedback on the back profile of the dune: they simply re- 
lax their equilibrium slope tant/j. Because no grains can 
escape the dune from the slip face through the recircula- 
tion bubble, the net out flux is zero, which fixes — Q. 
Note that this is particular to two-dimensional situations 
that we are focused on: three dimensional barchan dunes 
loose sand from their horns. 



6.1 The recirculation bubble 

The simplest way to model the feedback of the recircu- 
lation bubble on the whole dune has been proposed by 
Zeman and Jensen | [Tl| ] and used more recently by Kroy et 
al. Q . The idea is to build an envelope of the dune which 
prolongs the dune profile by the separation streamline - 
see figure |l^. To the first order, the wind on the back of 
the dune is the same as that would have been obtained if 
the envelope was solid. For example, the convolution inte- 
gral in equation (|^) used to calculate the shear stress on 
the soil should be applied not to the dune profile but to 
the dune -I- bubble envelope. 

In our approach, the effect of the recirculation bubble 
on the wind is simply to modify the total length of the 
dune D which becomes the dune length L plus the bubble 
length R: D = L + R. The importance of the recirculation 
bubble becomes very clear: because it makes the apparent 
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length of the dune be larger, it increases the erosion at the 
top of the dune, due to the curvature effect. 

We are going to construct a very simple description of 
the recirculation streamline. Because it extends smoothly 
the dune shape, it will be possible to link its parameters 
to the dune ones, and to impose this way what could be 
called 'non-local right boundary conditions' to the solution 
expressed by equation ([l9|). 

Let us call fefc(x) the separation streamline profile. As 
shown on figure O, x = L is the point at which this line 
starts and x = L + R is its reattachment point with the 
soil. At both matching points, the extension of the dune 
or the soil by the separation streamline should be smooth: 



hb{L) = h{L) = B, 
K{L)^h'{L)=p, 

hb{L + R) = 0, 

K{L + R) = 0. 



(30) 
(31) 
(32) 
(33) 



The recirculation bubble will be essentially governed 
by one important parameter: a slope that we already named 
Hb in the previous section. It corresponds to the maximum 
angle that the wind streamlines can follow behind an ob- 
stacle. In other words, there is a flow separation if the 
slope is locally steeper than —/if,. It is natural to impose 
that the steepest slope of the separation streamline should 
then be —^b- 



h'bixb) 



with 



:{xb) - 0. 



(34) 



As a matter of fact, this steepest slope should be in prin- 
ciple less than or equal to —fib otherwise this would mean 
that the streamlines actually detach at a higher thresh- 
old value. Besides, the bubble initially appears exactly at 
this slope on critical domes. The relation ( p4| ) can be then 
taken as the simplest self consistent slope condition. 

Because the wind velocity profile is scale invariant, a 
further property is that the length R of the recirculation 
bubble should scale on the dune size and no new length 
scale should be introduced by the bubble. Therefore, the 
curvature at x = L must also be continuous: 



h'^{L) ^ h"{L) = C. 



(35) 



Such a condition is also naturally observed when numer- 
ically integrating equations 
volution term of equation \\ 



(^,^,|5|) - i.e. with the con- 
J) rather than the simplified 
relation (|l^). Furthermore, we noticed that the stabil- 
ity of the numerical scheme is very sensitive to the way 
the slope is computed ai x — L: the only choice which 
does not create any instability is that which uses a for- 
mula which mixes the dune and bubble profiles h'{L) = 
[hb{L + Ax) — h{L)]/Ax, where Ax is the discretisation 
step. 

In the absence of available data or models on flow sep- 
aration behind a dune, the explicit form of the bubble 
profile we chose is the simplest that satisfies the require- 
ments written above, namely a polynomial of 3'^'* degree: 



1 



1 



B +p{x~ L) + -C{x- Lf + -G{x- Lf. (36) 
2 D 




-10 10 

horizontal position 



30 



Fig. 12. Dunes of different masses and corresponding recircu- 
lation bubbles. From small to large, the masses are M = 6.05, 
16.6 and 88.7. The other parameters are a = 1., /3 = 4. and 
/ifc = 0.25. For clarity, all these profiles have been shifted in 
order to get all brink positions at 1 = 0. 



Such a choice was also that of Kroy et al. in Using 
the conditions ( ^0[]35| ), G, C and R can be computed as 
functions of the slope p and the height B at the brink: 




2B 
3Rfib 

C = -RG + ^/2Gfib 
2G 



AB 
3Rfib 



P 



-l-i-b 



(37) 

(38) 
(39) 



We have checked that other similar choices for the parame- 
trization of the separation streamline do not change the 
qualitative conclusions that are going to be presented in 
the next subsections. In a general way, the recirculation 
bubble conditions can be expressed as R/B = fi{p) and 
CB — f2{p), where the functions /i and /2 encode the 
particular choice of the streamline separation profile hb. 

These two relations put together with the explicit ex- 
pression of the dune profile (19) and the differential equa- 
tion that it verifies ([l8|), where D has been set to the 
whole dune -I- bubble length D — L + R, let one to get 
three relations which link together the four parameters c, 
L, R and B. Instead of plotting, say, the three first with 
respect to the last one, we rather chose to express, as we 
did for the domes in the previous section, all of them as 
functions of the total mass M of the dune: 



M 



B 

dxh(x) + . 

2 tan (f 



(40) 



We then get a continuous set of dune solution, from very 
large values of the mass, down to some cut-off value Mc 
below which no stable recirculation bubble can be con- 
structed - see below. 



6.2 Profiles, dimensions and velocities 

Let us fix a value of the mass of sand M available to con- 
struct a dune. As for the domes, even for the simple choice 
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Fig. 13. Rescaled profiles of dunes of different sizes. The data 
have been computed with a = 1., (3 = 4. and fit = 0.25. 
Lengths have been rescaled by the square roots of the masses 
of the dunes. These profiles are not scale invariant: as shown 
in the inset, the slope p just before the brink is negative for 
small dunes and positive for large ones, such that depending 
on the dune size, the crest does or does not coincide with the 
brink. The four stars represent the four dune profiles. 




5 10 15 20 

dune size M"^ 



Fig. 14. Scaling of the lengths L, R (a) and H, B (b) with 
the dune size M^^^ . As a first approximation, these plots are 
almost straight lines. No recirculation bubble - and therefore 
no dune solution - can be constructed below some cut-off value 
(big dot). The inset on graph (b) is a zoom around this cut- 
off scale. Graph (c) shows that the dune and bubble aspect 
ratios are not constant: large dunes are more compact with a 
proportionally larger bubble than small ones. The data have 
been computed with a = 1., /3 = 4. and /ib = 0.25. 



3 




2 4 6 8 10 

dune height H 

Fig. 15. Propagation velocity c of the dune as a function of its 
height H. As shown in the inset, with a very good precision, 
1/c ~ aM^^^ + b down to some cut-off value. These data have 
been computed with a — 1., — 4. and nt = 0.25. 



of the recirculation bubble profile (pq), a numerical reso- 
lution is required at this point to get the corresponding 
values of the parameters c, L, R and B. Three examples 
of such solutions are plotted on figure As evidenced 
on figure ^ when rescaled by a typical dimension of the 
dune - here the square root of the mass -, the different 
profiles do not collapse on a single curve. In particular the 
slope at the brink p = h'{L) varies with respect to M, and 
even changes its sign for this choice of the parameters a, 
(3 and fib. We get p < for the smallest dune, while p > 
for the largest one. 

The inset of figure as well as the curves of figures 
|l^ and |l5|, which all show the dune and bubble features 
(lengths, slope, aspect ratios and velocity) as a function of 
the mass of the dune M, are all cut off at some value Mc 
(big dot). This point corresponds to the smallest solution 
that can be constructed. Technically, one can see that the 
equations ( ^7| - |39| ) which link bubble and dune parameters 
loose their sense when the quantity SRfib — 2-B becomes 
negative. This point precisely gives the cut-off mass un- 
der which no bubble, and therefore no dune solution can 
be constructed. It is important to notice that this small- 
est dune has an avalanche slip face of finite size. Another 
choice for the function hb{x) would of course have given 
different values for Mc as a function of the parameters of 
the model. But physically, this cut-off scale corresponds 
to the fact that such a recirculation bubble must have a 
minimum spatial extension to accommodate all continuity 
constrains. Not surprisingly, this extension is of the order 
of unity, that is to say of the order of Igat in non rescaled 
length units ~ which is the only length scale of our system. 
At last, it must be mentioned that the parameters of this 
smallest dune are pretty close to that of the largest dome. 
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M 


c 


H 


Smallest dune 


4.51 


2.78 


0.86 


Largest dome 


6.55 


2.78 


0.78 



Table 1. Mass, velocity and height of the smallest dune so- 
lution that one can construct, compared to that of the largest 
dome for which the steepest slope is equal to — fib - Consistently, 
these data are pretty close to each other. These numerical val- 
ues have been computed with a = 1., /3 = 4. and Hb = 0.25. 



0.2 



0.25 



parameter 
0.3 



0.35 



— H>,; 



parameter ii,, 
0.3 



£ 0.1 



0.4 g_ 




parameter (3 

Fig. 17. Variation of the asymptotic (M —f oo) slope at the 
brink p with the parameters /3 and fib- To plot theses curves, 
we took a = 1- Besides, the solid line has been obtained with 
fib = 0.25, and the dashed one with P — 4. Interestingly, p 
remains negative for small /3 or large fib- 



parameter p 

Fig. 16. Variation of the minimal dune height He with the 
parameters /3 and fib- To plot theses curves, we took a = 1- 
Besides, the solid line has been obtained with fib = 0.25, and 
the dashed one with /3 = 4. 



i.e. the dome for which the steepest slope is precisely equal 
to —fib, see table |l|. 

All results of figures and 15 are remarkably simple 
and resemble very much field observations: the lengths L, 
R, H and B are almost straight lines as a function of the 
dune size M^^^. Similarly, 1/c ~ aM^/^ -I- b with a very 
good precision. At last, it must be noted that the dune 
and bubble aspect ratios are not constant: large dunes are 
more compact with a proportionally larger bubble than 
small ones. 



6.3 Parametric study 

To complete the results of the previous subsection, we 
present here a parametric study of few quantities, namely 
the minimal dune height He and the dune slope at the 
brink of asymptotically large dunes {M — > oo). The pa- 
rameter a is kept to unity, but /? has been varied between 
2 and 6, and fib from 0.2 to 0.4. 

Figure |6| shows that the variation of He as a func- 
tion of the parameters of the model is not very strong. 



Similarly, the aspect ratio of very large dunes is always of 
the same order of magnitude. The most interesting fact is 
perhaps that for small values of /3, or for large values of 
fib, the slope at the brink of very large dunes can remain 
negative - see figure 17| In other words, the brink and the 
crest can be always distinct, even for asymptotically large 
dunes. 

One can understand intuitively the variations of He 
and p with /3 and fib- Increasing gives more strength to 
the destabilizing process, which lets small dunes appear 
at lower critical scale (smaller He) and makes large dunes 
more bumpy (larger p). If fib gets larger, the situation is 
reversed. 



7 Conclusion 

We have shown in this paper how, inspired from the work 
of Sauermann et aL, one can build a simpler two dimen- 
sional model for the formation and the propagation of 
dunes. This modelling is based on two main variables: the 
dune profile h and the volumic sand flux q. It includes 
three effects: (i) the mass conservation, (ii) the space lag 
over a length Isat for the sand flux to become saturated 
at some value qsat, and (iii) the feedback of the profile on 
the saturated fiux. In this third phenomenological equa- 
tion, erosion and deposition processes are the result of the 
competition between two antagonist mechanisms: a sta- 
bilizing non local curvature term (a) and a destabilizing 
slope one (/3) which breaks the upwind-downwind symme- 
try. 
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Two kinds of solutions have been found: so-called 'domes 
which do not show any avalanche slip face, and 'real' dunes 
for which a downwind recirculation bubble has been in- 
troduced. We were able to predict an analytical form for 
their propagative profiles, but whose coefficients have to 
be computed numerically. 

The results of the model resemble very much field ob- 
servations. We found for example that, due to this length 
scale I sat, the dune profiles are not scale invariant: small 
dunes are fiatter than large ones. Another point is that 
the inverse of the propagative velocity c is, to a very good 
precision, almost linear with the size of the dune. This is 
consistent with Bagnold's argument that c ~ Qsat/H for a 
dune of height H. In fact, this relation overestimates the 
velocity of small dunes for which the sand flux may not 
be already saturated at the crest and also because, com- 
pared to large dunes, the value of qsat is reduced due to a 
smaller curvature. 

An important point discussed all through the paper 
was the issue of the boundary conditions. In particular, 
an important physical input was the 'recirculation bubble' 
behind the dune. This bubble makes the dune effectively 
look larger to the wind, and, due to the non local term 
in the relation between the saturated fiux and the dune 
profile, has a stabilizing role. In fact, very little is known 
and well established about this recirculation bubble, but 
most of the dune features (position of the slip face, cut- 
off size, etc.) precisely depend on fine interactions and 
feedbacks between the dune and the bubble. More studies 
on this point are then urgently needed. 

Another central result of the paper is the existence of 
this cut-off scale just mentioned, below which no dune can 
form. It corresponds to the fact that the bubble must have 
a minimum spatial extension to accommodate all conti- 
nuity constrains. Not surprisingly, this scale is of order 
of I sat- This result then rises the question of dune initia- 
tion and formation. The two scenari usually proposed by 
geophysicists for the formation of dunes are the following: 
first possibility, a small bump (of the size of ripples) grows 
continuously and forms a dune; second one, the sand ac- 
cumulates on a solid obstacle like a rock or a bush and, 
when the size of the accumulation becomes larger than 
the obstacle, a dune forms and starts propagating down- 
wind. However, observations show that ripples are stable 
and no structures between dunes and ripples can be seen. 
Similarly, rocks and bushes creates lee dunes of the size 
of the obstacle which remain anchored to the obstacle. 
An alternative explanation can be proposed, following the 
results of the stability analysis of the equations of the 
model, as well as that of the dome solutions. We found 
that large wavelengths perturbations get amplified, and 
that the dome profiles, selected by their incident flux qq, 
are unstable to changes of that flux. Then, a possibility 
is that first domes form with a small height but directly 
with large length and width, and second that these domes 
progressively become more and more compact, and even- 
tually reach the point where their slope is steep enough 
to generate a bubble and create an avalanche slip face to 
become an actual dune. 



Several extensions to the present work can be thought 
of. First, we would like to go beyond the calculation of 
purely propagative solutions, and study the full dynam- 
ics of a given dune profile. In particular, as just said, an 
important point is the evolution of the dome solutions 
when submitted to incident sand flux variations. A second 
point is to go from a 2d description - transverse dunes - 
to real three dimensional situations. The idea is to 'cut a 
barchan into longitudinal 2d slices'. As a matter of fact, 
a barchan slice close to the center of the dune looks like 
our dune solution, while a slice made at the edges where 
the horns are present rather have a dome shape. Suppose 
these slices are completely decoupled. Because the small 
ones go faster than the large ones, an initial conical sand- 
pile will soon get a crescentic shape. However, when equi- 
librium is reached, all the slices should move at the same 
velocity. There should thus be a coupling between them, 
namely a lateral sand flux from the centre towards the 
horns. When the flux is saturated at the crest the veloc- 
ity at the crest is c = {qsat — qo)/H. This suggests than 
an equilibrium can indeed be achieved if go increases in 
the small slices. Eventually the 3D dune slip face will be 
the sum of the contributions of all 2D slices whose brinks 
depend on their heights. Note that this scenario is con- 
sistent with the field observation that barchan horns are 
more elongated at strong winds which make the lateral 
sand flux less important, and consequently slices less cou- 
pled. 

Finally, quantitative comparisons between experimen- 
tal dune proflles and our theoretical predictions will be 
performed. This idea is to use barchan longitudinal slices 
as that shown on flgure |^, but also sand structures under 
water, such as those obtained by Betat et al. or Ander- 
sen et al. [p^ for which, in principle, this model should be 
also valid. 
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